NUMERICAL ASPECTS OF NONLINEAR SCHRODINGER 
EQUATIONS IN THE PRESENCE OF CAUSTICS 



REMI CARLES AND LAURENT GOSSE 

Abstract. The aim of this text is to develop on the asymptotics of some 1-D 
nonlinear Schrodinger equations from both the theoretical and the numerical 
perspectives, when a caustic is formed. We review rigorous results in the field 
and give some heuristics in cases where justification is still needed. The scat- 
tering operator theory is recalled. Numerical experiments are carried out on 
the focus point singularity for which several results have been proven rigor- 
ously. Furthermore, the scattering operator is numerically studied. Finally, 
experiments on the cusp caustic are displayed, and similarities with the focus 
point are discussed. 



1. Introduction 

We present a numerical study of the semi-classical solutions to the following non- 
linear Schrodinger equations with e <C 1, 

(1.1) ieatu" + yAu^ |u'^P'^u% (t,a;)eM+xM" ; uf^^g = eP/(a;)e*"^''(^)/% 

when a caustic (a point or a cusp) is formed, that is to say, beyond breakup time. 
Since the nonlinearity is homogeneous, the change of unknown function u'^ = e^u*^ 
shows that (jl.ip is equivalent to: 

(1.2) ledtu' + —Au' ^e^^Plu'l^-^u' ; uf^^o = /(x)e*■^«(^)/^ 

so that we can always consider initial data of order 0(1). 

There are several motivations to study the behavior of p.2p when a caustic is 
formed. First, on a purely academic level, we recall that the description of the 
caustic crossing is complete in the case of linear equations; see |17j . For nonlinear 
equations, very interesting formal computations were proposed in [25j (we recall 
the main idea in Section [2] below) . For dissipative nonlinear wave equations, Joly, 
Metivier and Ranch [551 127] have proved that the amplification of the wave near the 
caustic can ignite the dissipation phenomenon in such a way that the oscillations 
(that carry highest energy) are absorbed. The above nonlinear Schrodinger equation 
is the simplest model of a conservative, nonlinear equation. The mass and the 
energy of the solution are independent of time (see (12. 7|) below). Therefore, different 
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nonlinear mechanisms are expected. We recall in Section [2] some results that have 
been established rigorously, and give heuristic arguments to extend these results. 
This serves as a guideline for the numerical experiments proposed after. 

Second, (|1.2p may be considered as a simplified model for Bose-Einstein conden- 
sation, which may be modeled (see e.g. [I6l[32|) by: 

(1.3) iedtu' + yAw" = + 

with (T = 2 if n = 1, and cr = 1 if n = 2 or 3. The power in front of the 
nonlinearity depends on the regime considered, and in particular on the respective 
scales of different parameters (see e.g. [7] and references therein). The role of 
the harmonic potential \x\'^ is to model a magnetic trap. In the semi-classical 
limit e — > for the linear equation, this potential causes focusing at the origin for 
solutions whose data are independent of s. This is to be compared with the case 
of (jl.2p with initial quadratic oscillations as considered below: the initial quadratic 
oscillations force the solution to concentrate at one point in the limit e ^ 0. The 
parallel between (II. 2p and (II. 3|) was extended and justified in [10] for these nonlinear 
equations. 

From both points of view, when a caustic point is formed, the caustic crossing 
may be described in terms of the scattering operator associated to 

This aspect is recalled in Section [21 For this reason, we also pay a particular 
attention to this operator, independently of the above semi-classical limit. Note 
that besides the existence of this operator, very few of its properties (dynamical, 
for instance) are known. 

In this paper, we always assume 2ap > 1: one of the reasons is that when 
< 2ap < 1, instabihty occurs, see [HI [13]. Suppose for instance that solves 
(|1.2p . and that solves (jl.2|) . where / replaced by {1 + S'^)f, where S'^ is a sequence 
of real numbers going to zero as e — > 0. Then there are some choices of for which 

limM\\u'{f) -u'{f )\\L2 > 0, 

for some sequence of time ^ (see flT , and f7] for a similar phenomenon with 
different initial data). Therefore, producing reliable numerical tests in the case 
< 2ap < 1 (which is super-critical as far as WKB analysis is concerned [l2l [l3] l 
seems to be a very delicate issue, that we leave out in the present paper. 

The rest of this paper is structured as follows. In Section [51 we recall the general 
approach of WKB analysis for the Schrodinger equation, the arguments of [25] . 
and the rigorous results available for the semi-classical limit of (|1.2p when a caustic 
reduced to a point is formed. We then recall the definition of the scattering operator. 
We also give heuristic arguments to tackle the case of a "supercritical focal point" , 
and to guess what the critical indices are when a cusp caustic is formed, instead 
of a focal point. In Section [31 we present the different strategies that have been 
followed in the literature to study numerically the semi-classical limit for nonlinear 
Schrodinger equations. Numerical experiments on the semi-classical limit for (|1.2p 
in the presence of a focal point appear in Section [H and the scattering operator is 
simulated in Section [5] We present the numerical experiments of the semi-classical 
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limit for ()1.2|1 in the presence of a cusp caustic in Section [6l and make conclusive 
remarks in Section [T] 

2. Analytical approach 

2.1. Semi-classical limit of the free Schrodinger equation. Consider the ini- 
tial value problem, for {t,x) £ R+ x M": 

(2.1) isdtv' + '-^Av' ^0 ; z;f,^o = /(a;)e'*''(-)/^ 

The aim of WKB methods is to describe the asymptotic behavior of w*^ as £ — > 0. For 
instance, e can be related to the Planck constant, and the asymptotic behavior of 
is expected to yield a good description of when e is fixed, but small compared 
to the other parameters. More precisely, seek of the form 

(2.2) v'{t,x)r^e"t'^'^'''>^''{aQ{t,x)+eai{t,x) + ...) as £ ^ 0. 

Plugging this expansion into (|2.ip and canceling the 0{e''^) term, we see that the 
phase 4> niust solve the eikonal equation: 

(2.3) 5,0+i|V(/.p=O ; (l)\t=o = ^o- 

To cancel the 0{e^) term, the leading order amplitude solves the transport equation: 

(2.4) atflo + V(/) ■ Vao + iaoA0 ; ao|t=o - /. 

The eikonal equation (|2.3p is solved thanks to Hamilton- Jacobi theorjQ: 4> is con- 
structed locally in space and time (see e.g. ^3j for a discussion on this aspect). 
Even if (po is smooth, (j) develops singularities in finite time in general: the locus 
where </> is singular is called caustic (see e.g. the second volume of 124j). When (f> 
becomes singular, all the terms ao,ai, . . . may become singular as well. One easily 
observes that (j2.4p admits a "divergence form": 9t|aop + V • (|aopV0) — 0. To 
illustrate this general discussion, we consider two examples that will organize the 
rest of this paper. 

I |2 

Example (Quadratic phase). Let (t>o{x) = Then p.3p and (|2.4p can be 

solved explicitly: 

0(i x) = -J^i— ; a(t,x) = - ^-^^f(^^\ 

This shows that a.s t 1, (j> and a become singular: the wave focuses at the 
origin. This example can be viewed as the smooth counterpart of the Cauchy 
problem 

1 

idti^ + -^A^j = ; V|t=o = e ' = . 
Fourier analysis shows that "014=1 = 6, the Dirac measure at the origin. 

Of course, the solution of (|2.ip can be represented as an oscillatory integral: 

1 / -la; — vl^ , ■ din (y) 

(2-5) -'(t^-) - (^-^ J e^^^~f{y)dy. 



but not according to the theory of viscosity solutions! See e.g. | 28| . 
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The caustic set is exactly the locus where the critical points for the phase 

*M(2/) = ^^ + '^o(y) 

are degenerate. Outside the caustic, an approximation of is given by the sta- 
tionary phase theorem (that we recalled as simply as possible in [22]). This leads 
us to the second example we shall consider numerically: 

Example (Cusp). Let n = 1 and (t>o{x) — cos a;. The set of degenerate critical 
points for ^t.xiu) (caustic) is given implicitly by: 

[ u — X 1 

C = <^ (t, x) e K+ X K; Ely e M, ^ — = siny, and - = cos y 

As soon as t > 1, a caustic is formed (see Figure 2 in [53]). 

When considering the asymptotic behavior of beyond the caustic, two main 
features must be considered: the creation of other phased, and the Maslov index 
(see [T7] for more general linear equations). In the case of a focal point, the first 
aspect does not exist: there is no creation of phase, and one phase is enough to 
describe past the focal point {t,x) = (1,0). One can prove easily the following 
result: 

Lemma 2.1. Letn > 1 and f £ iS(]R";C). If(po{x) — —\x\'^/2, then the asymptotic 
behavior (in L^(M")j of the solution to (|2.ip is given by: 

j|x|V(2£(t-l)) / X \ 



^ - -0 1 ^ i|a;|V(2e(t-l)) 



(t-l)"/2 •' Vl-t, 

In this example, the Maslov index is — n7r/2. In the case of the cusp, three phases 
must be considered to describe the asymptotic behavior of beyond the caustic 
(see e.g. ^^). 

For future discussion on the numerical results, we state the following more precise 
result, which follows from the stationary phase theorem: 

Lemma 2.2. Letn > 1 and f € 5(]R";C). If 4>q{x) — ~\x\'^/2, then the asymptotic 
behavior of the solution to (|2.ip at time t = 2 is given by: 

v%2,x) = e-^fe^l^l'/^^'^)/ (-x) + 0(e) in n L°°(R"). 

2.2. Caustics in the nonlinear case: heuristics. Consider now the perturba- 
tion of (12. ip with a nonlinear term: 

(2.6) ledtu' + —Au' ^e^'lu'l^'^u' ; u^^^^ ^ f{x)e"t">^''^/' . 

The sign of the nonlinearity is chosen so that no finite time blow-up occurs. The 
following two important quantities are formally independent of time: 

Mass: \\u^{t)\\L2 = Const. = ||/||l2. 
^^■^^ Energy: E%t) := l\\eVu^t)\\l. + ^\\u^t)\\l%t^, = E^O). 

^in our mind, phases are always associated to oscillations whose period depends on e, and goes 
to infinity as £ — » - rapid oscillations. The wavelength may be proportional to e, or, say, to y/e. 
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We refer to [M] for a justification. Fix the power 2(7 > of the nonhnearity, 
and consider different values for a. Two notions of criticality arise: for the WKB 
methods on the one hand, and for the caustic crossing on the other hand. This 
discussion is presented in [25^ for conservation laws, and we summarize it in the 
case of (|2.6p . Plugging an expansion of the form (|2.2p into (|2.6p . we see that the 
value a = 1 is critical for the WKB methods: if a > 1, then the nonhnearity does 
not affect the transport equation (|2.4p ("linear propagation"), while if a = 1, then 
the nonhnearity appears in the right hand side of (|2.4p ( "nonlinear propagation" ) . 
Recall that in this paper, we always assume a > 1. Therefore, the eikonal equation 
(j2.3p is not altered: the geometry of the propagation remains the same as in the 
linear WKB approach, and we have to face the same caustic sets. The idea presented 
in [25] consists in saying that according to the geometry of the caustic, different 
notions of criticality exist, as far as a is concerned, near the caustic. In the linear 
setting (|2.ip , the influence of the caustic is relevant only in a neighborhood of this 
set (essentially, in a boundary layer whose size depends on e and the geometry of 
C). View the nonhnearity in (|2.6p as a potential, and assume that the nonlinear 
effects are negligible near the caustic: then u'^ ~ near C. View the term e"|u'^p°' 
as a (nonlinear) potential. The average nonlinear effect near C is expected to be: 



where is the region where caustic effects are relevant, and the factor e^^ is due 
to the integration in time (recall that there is an e in front of the time derivative 
in (|2.6p ). The idea of this heuristic argument is that when the nonlinear effects are 
negligible near C (in the sense that the uniform norm of — is small compared 
to that of v"^ near C^), the above approximation should be valid. On the other 
hand, it is expected that it ceases to be valid precisely when nonlinear effects can 
no longer be neglected near the caustic: — is of the same order of magnitude 
as in L°°{C^), or even larger. 

Practically, assume that in the linear case, v'^ has an amplitude e^^ in a boundary 
layer of size e'^; then the above quantity is 



The value a is then critical when the above cumulated effects are not negligible: 

etc = 1 + 2^0- - k. 

When a > ttc, the nonlinear effects are expected to be negligible near the caustic: 
resuming the terminology of |25j , we speak of "linear caustic" . The case a = ac is 
called "nonlinear caustic" . To conclude this paragraph, we examine this approach 
in the case of our two examples. In the case of a focal point, we have k = I and 
i = n/2. This leads us to the value: 



In the case of the cusp in dimension one, we have fc = 2/3 and = 1/3 (which 
can be viewed thanks to the Airy function and its asymptotic expansion, see e.g. 
[m EH [25] or [29]), which yields: 





ttc (focal point) 



= na. 



ac(cusp in ID) 



2(7+1 

3 



6 



R. CARLES AND L. GOSSE 



One aspect of the numerical experiments presented below is to test this notion of 
criticality in those two examples. 

2.3. Justification for a focal point, and more heuristics. In this paragraph, 
we assume 4)q{x) = —\x\'^/2. A complete justification of the above discussion is 
available [8]: 





a > na 


a = na 


a > 1 


Linear caustic, 
linear propagation 


Nonlinear caustic, 
linear propagation 


a = 1 


Linear caustic, 
nonlinear propagation 


Nonlinear caustic, 
nonlinear propagation 



Consider t E [0, 2], which includes the caustic crossing. The above tables means: 

• If a > max(l, na), then u*^ can be approximated by for t G [0, 2]. 

• If a = 1 > na, then the nonlinearity is negligible near the focal point, but 
not away from it. 

• If a = na > 1, then nonlinear effects are relevant near the focal point, and 
only near the focal point. 

• If a = na = 1, then the nonlinearity is never negligible. 

We give some precisions in some cases of interest for the numerics presented below. 



In the one-dimensional case n — 1, the following pointwise estimate is proved in 
[8| when a > max(l, ct) or a = ct > 1: 

C 



(2.8) 
Setting 



— f ^ , we see that 

^2 



The usual energy estimate yields, for t > 0: 



\\w%t)\\L2 < e-' / e'^\\\u%r)\^'^u%T)\\^,dT. 
Jo 



Using ([23]), we infer, for t G [0, 2]: 



dr 



{|r-l|>e}n{r6[0,2]} F " J\T-l\<e 



dr 



Using the operator edx and x/e + i{t — l)dx, and Gagliardo-Nirenberg inequalities 
as in [8], we find: 

Lemma 2.3. Let n = 1, a > max(l, ct), / G 5(R"; C), and (?!)o(a;) = -|a;p/2. Then 
we have, for the solutions of (j2.ip and (|2.6p ; 



sup \\u''{t)^v'{t)\\ < Ce"-" 



\u%t)~v''{t)\ 



Q<t<2 

In particular, Lemma \2.3\ implies, at time t — 2: 



^\t-l\+e 



, te[0,2]. 
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(The above result is true also when a = cr > 1, but becomes far less interesting.) 
We now explain the critical case a — na > 1. The nonlinear effects near the focal 
point are described in terms of the scattering operator associated to the nonlinear 
Schrodinger equation. We rapidly present this operator S in Section \ZM We have 
then: 

{t - i)"/2 VI - V 

where Z = o S o T^^ is the conjugate of 5* by the Fourier transform (see [5]), 

(Since S* is a nonlinear operator, the normalization of the Fourier transform is an 
important detail.) 

To conclude this paragraph, we give a few hints of what happens or may happen 
when the propagation is linear, and the caustic is super-critical, that is 1 < a < na. 
First, the conservation of mass and energy seem to rule out the possibility of a 
concentration of the form 

(2-9) -^^(l'-)-^K?)' 

for some function Lp independent of e. The above relation holds for and when 
a > max(l, ncr), with kp — T j (for w^, this is obvious from (j2.5p ). When a — na > 1 
(linear propagation and nonlinear focal point), the above relation still holds, with 
a different profile ip (see |8]). Now we see that the energy is bounded as e — > 0: 

Plugging a concentrating profile as in (|2.9p in the second term of the energy would 
yield, thanks to the conservation of mass: 

e"h^(l)||^'^+^«e™— 

which is impossible since the energy is the sum of two positive terms. This suggests 
two possible effects: near t — 1, the amplification of the solution (in terms of 
powers of e) may be weaker than in the linear case; on the other hand, nonlinear 
effects near the caustic should affect the phase of the solution in a rather strong 
way, causing the appearance of new frequencies. A partial justification of the last 
assertion may be found in [11] . 

Extrapolating this argument, we expect that in the supercritical case for a cusp 
(with n — 1: > a > 1), new frequencies appear. If as in the linear case, three 

phases are necessary to describe the solution past the caustic, then the nonlinear 
interaction of these phases might reveal the presence of new frequencies, even on the 
modulus of M*^ (see Sect. H] for numerical tests that seem to confirm this heuristics). 



8 



R. CARLES AND L. GOSSE 



2.4. The scattering operator for NLS. To explain what the operator S men- 
tioned in the previous section is, consider the nonhnear Schrodinger equation 

(2.10) idti^ + ^ Ail; =1^1^" ij ; {t,x)eRxW\ 

(2.11) U{^t)m\t=to=^~, 

where U{t) — e'^'^ is the propagator of the hnear equation. To construct the 
scattering operator, we first want to give a meaning to (|2.10p - (l2.1ip when to = ~oo. 
This means that the nonhnear effects are asymptoticahy negligible as t — > — oo: for 
instance, we expect at least 

\\u{-t)ij{t) - v_ = wm - iil^ 0. 

t — * — oo 

This gives a rigorous meaning to the relation ip(t,x) ~ U{t)'ip-{x) which aims at 
saying that as time goes to — oo, the nonlinear dynamics associated to (|2.10p can 
be compared to the free dynamics given by e'^^. 

To define a scattering operator, we want to be able to say that as t — > +00 
as well, the nonlinear effect are asymptotically negligible. That is, there exists 
4'+ e L2(M") such that 

\\ui^t)m - V'+iu^ = wm - (7(ov+iu2 — . o. 

t — ^+00 

The scattering operator is then defined as S : '0_ >—>■ ?/'+ • Since our numerical exper- 
iments concern the one-dinicnsional case, we recall the existence of the scattering 
operator in this setting, and refer to [T31 [HI [201 [HI [31] for some extensions to the 
multidimensional framework: define 

S := H'nTiH') = If e L2(M") ; ll/lls := + + ||V/|U. < 00}. 

Proposition 2.4 (Scattering theory). Let n = 1, and assume a > ^^^^ {> 1). 

• For every tjj- £ S, there exists a unique 95 G S such that the maximal 
solution ijj G C(R, S) to ()2.10p satisfies V'|t=o = f o.'^d 

\\u{-t)m-mk,^ 0. 
t— >— 00 

• For every £ E, there exists a unique G S such that the maximal 
solution ijj G C(M, S) to p.lOp with V'|t=o — V satisfies 

||[/(-i)^(t)-V^+|l5, — > 0- 

The scattering operator is S : ip- t— > ip^. When a > 1, the above conclusions 
remain, in a neighborhood of the origin. When a > 1, we also have: 

• For every G S, there exist a unique solution ip G C{M.,H^) to (|2.10p 
and a unique ip+ G H^{M) such that: 

\\u{-t)m-i^-\k,^ ^ ; |jc/(-t)VXi)-V'+llL2 — ^ 0. 

When a < 1, the above conclusions are false: if cr = 1 for instance, and if 
tA- G with Uo{-t)ii{t) --0-^0 in as t ^ -00, then ^ = 0. One 

cannot compare the nonlinear dynamics with the free dynamics (see [3l [33ll34lll9) V 
Note that even though the scattering is proven to exist, very few of its features are 
known. We refer to for some algebraic properties. At least, this operator is not 
trivial: near the origin, it is a non-trivial perturbation of the identity (see [5]). 
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3. Numerical approximation of semi-classical Schrodinger equations 

Hereafter we restrict our discussion to the one-dimensional case, that is to say 
n = 1 in all the preceding considerations. 

3.1. Rigorous results for general time-splitting schemes. It is interesting 
to notice that in the numerical literature, the nonlinear equation (j2.6p is treated 
exactly the same way the linear one (j2.ip would be in the presence of a potential 
V on its right-hand side. The strategy is called time-splitting, in its first or second 
order version (Lie or Strang splitting, see e.g. [1]) where one alternates every time 
step At > between the solving of the Laplace operator and the handling of the 
(nonlinear) differential equation. According to Section 12.41 U will still stand for 
the free propagator, whereas we shall use V as the ODE solver; Lie time-splitting 
algorithms generate the following type of approximation for (|2.6p . 

u'inAt, .) ~ u%^{nAt, .) := [V{At) o f/(At)]"u"(i = 0, .), 

and Strang splittings, 

u^{nAt, .) ~ UAtinAt, .) V{At/2)oU{At)[V{At)oU{At)]''-^V{At/2)u^{t = 0, .). 

Many references exist; let us quote only [4[ l30l [Tl [2]. 

On the contrary, few rigorous convergence results are available, hence we shall 
mainly recall the results from [U which quantify accurately the splitting errors 
assuming each time-step is performed exactljQ. Under this assumption, there holds: 

Proposition 3.1. ([4^, Theorem 4-1) For any T > and (t 0, .) e iJ^, there 
exists a constant C depending on the initial data for (j2.6p and Hq, such that for 
At G [0, ho] and uAt < T, 

\\u'{nAt, .) - u'^tinAt, .)\\l2 < Cho- 

If moreover u^{t = 0, .) G , then there holds under the same assumptions: 

\\u'{nAt,.)-u%t{nAt,.)\\L2<Chl. 

This result is concerned with splitting errors only and relies on the knowledge 
of the exact solution operators U and V . In order to stick to this framework in 
the context of smooth solutions, it is rather natural to approximate U by means of 
a Fourier scheme taking advantage of optimized FFT routines, as proposed in the 
paper f2j. Moreover, this will guarantee that the norm (Mass) of the numerical 
solution will be conserved up to round-off errors. Unfortunately, the Hamiltonian 
E'^it) is generally not preserved; a method conserving both quantities exists (see 
the so-called MCN algorithm, page 253 of [18]) but it wouldn't be efficient in the 
semiclassical regime because of the results in |30j . 

3.2. Specific issues witli finite-difference discretizations. This is the main 
purpose of the paper [SOj to illustrate the (surprising) fact that in semi-classical 
regime, usual finite-difference schemes for (j2.ip can deliver very wrong approxi- 
mations without any particular sign of instability in case very restrictive meshing 
constraints turn out to be bypassed. This can be quite easily checked through the 
location of caustics, for instance. The analysis of those standard schemes has been 
carried out by means of Wigner measures, so the conclusions hold essentially for 
the quadratic observables coming out of the wave function itself. 



^but we shall see in the sequel that this is far from being the case! 
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3.3. The case of FFT-based schemes. This class of schemes became popular 
after the publication of [HE] , mainly because treating the differential part of (12. 6p by 
means of a discrete Fourier transform looked very much like being the best possible 
compromise in terms of meshing constraints. Indeed, in the linear case where (12. 1[) 
is supplemented with a potential V{x) on its right-hand side, it was shown that 
the time-step At could be chosen independent of e whereas the space discretization 
has to satisfy Ax = 0(e). This was already much better when compared to finite- 
differences; moreover, the method is naturally L^-conservative. In [2], these authors 
extended their "Fourier framework" to the weakly nonlinear Schrodinger equations 
of the form (|2.6p . However, and despite the fact we do believe these "FFT time-split 
schemes" realize the best numerical strategy in terms of gridding, we shall point 
out some shortcomings of the method in the next section. 

3.4. A rigorous framework for FFT-based schemes. We present here a pre- 
liminary result about truncation errors in Lebesgue spaces for Fourier schemes; its 
proof follows directly from the Strichartz estimates on the torus due to J. Bour- 
gain [5] (see also 0), and from the study of FFT by M. Taylor [35l pp. 250-254]. 
Its derivation is not obvious though as it applies directly to widely-used schemes 
like the one recalled in the forthcoming section. We restrict our attention to the 
ID free Schrodinger equation (|2.1|) with e = 1, and periodic boundary conditions: 
xGT := R/27rZ. 

Hence we start from 

19*^ + ^52^ = 0, V(i = 0,-) = C = EOe*'". xe[0,27T]. 
We have explicitly: 

V'(t,x)=5]0e^^(--^"*/2). 

In order to investigate the behavior of the FFT-schcmc involving a finite even num- 
ber N G 2N of modes, we introduce the Discrete Fourier Transform of a continuous 
function / on [0, 27r] as follows: 

From 35, p. 252], we recall: 

Lemma 3.2. // the continuous function f has a convergent Fourier series, then: 

Now, what the Fourier numerical scheme really computes is: 
Thus the corresponding numerical solution ^^^"^"^ is built: 

N_ 
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Num. . 



The question is therefore to study the discrepancy between ip and ip 



jez feez 



:(a;-fei/2)A 1 ik(x-kt/2) 





^( 


E 


jez \ 




E 

jez* 





A ik(x-kt/2) 
,ke ^ 



|fc|>- 



The second term // in this last equahty can be interpreted as a usual truncation 
error; it means that some high frequencies in the solution are lost when discretizing 
the problem on a finite grid. From |351 p. 253], we know that the modulus of this 
term can be controlled as follows: 

E iai<^iiciic^+^, 
ifci>f 

which is satisfactory provided the initial data C is smooth, e.g. C^(R). However, 
the first term / is far more delicate and reveals the propagation of the errors coming 
from the DFT/FFT itself. It can be controlled thanks to the periodic Strichartz 
estimates proved by J. Bourgain. From [6, pp. 16-17], we recall: 

Lemma 3.3. Given any complex-valued sequence {ak)k, the following estimates 
hold for the corresponding functions in [0, 27r] x [0,27r]; 

||E«fce'^''"''*^ ' <cEk^ 

II ^ L^fTxT) ^ 

k k 

\k\<N ^ ' k 

The second estimate looks more attractive as it "sees" the finite number of 
modes. We therefore deduce that the first term can be controlled in LP' by means 
of: 

^ ^ ^ I ^ VjezMfc|<f 

For instance, if is a finite superposition of Fourier modes, then it is clear that this 
term cancels for N large enough as 7^ in the summation; obviously, the second 
term // vanishes too. 

The general case isn't completely clear yet. 

4. Experiments on the focal point 

This section aims at visualizing the asymptotics previously recalled; namely we 
shall compare numerical approximations of (|2.6p and (|2.ip in ID (n = 1) for various 
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values of the parameters a and a. The initial wave function is rather simple: 

v'it = 0,x)= u''{t = 0,x) = exp (^-(2 + -tt)^^ , x e [0,2n]. 

Numerical results have been obtained through the time-splitting FFT schemes re- 
called in the previous section; we used 1024 modes and fixed e = 1/150. It is 
convenient to observe results in t = 2 since \v^{t = 2, .)| = \v'^{t — 0, .)|. 

4.1. Subcritical case. This case corresponds to cr = 2 and a = 2.5; we expect to 
observe a decay of the absolute errors between v'^ and for values e <C 1. This is 
indeed the case, but Fig. [T] shows even a bit more, namely it compares pointwise 
the following quantities: (recall Lemma [2T|) 

SR (v%t = 0, x) exp(-i(a; - T:)^/2e)) , 
S [v^lt = 2, x) exp(i(x - tt) V2e)) , 
S (u^(t = 2, x) exp(i(a; - Tr)^/2e)) . 

On the left in Fig. [U we obviously observe that the absolute errors are slightly 




Figure 1. Absolute errors on the wave functions (left) and on the 
modulus (right) at T = 2 with a = 2 and a = 2.5. (Subcritical 
case) 

bigger when considering the solution of the nonlinear equation (j2.6p . . However, 
even for the free solution u^, one sees that the error doesn't vanish despite the 
fact no time-splitting algorithm is needed. As the way of discretizing the solution 
reveals itself important, we include here the corresponding Scilab routine for the 
free equation: 

clear ; def f ( ' [y] =phase (x) ' , [ ' y=-0 . 5* (x-tt) ' ; ] ) 

deff (' [y]=position(x) ' , ['y=exp(-2*(x-7r) .2) '] ) 

def f ( ' [y] =Az(x) ' , ['y=position(x) . *exp(i*phase (x) . / epsilon) ' ] ) 

NMAX=2i" ; n=- (NMAX) /2 : (NMAX/2) -1 ; epsilon=l . 0/150 ; 

XSTART=0 ; XSTDP=2*7r ; DX= (XSTDP-XSTART) /NMAX ; XSTOP=XSTOP-DX ; 

a=XSTART : DX : XSTDP , initialdata=Az (a) ; 

vepsilon=f f tshif t (f f t (initialdata, -1) ) ; 

vepsilon=exp(-i*epsilon*(n.^) ) . *vepsilon; 

vepsilon=f f t (ff tshif t (vepsilon) , 1) ; 
Clearly, its outcome is in agreement with Lemma |2 . 1 1 since e is already quite small. 
The Maslov index is visible, up to an error around 10^"^ for 1024 Fourier modes. 
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4.2. Critical case. We now put a ~ a = 2 and the outcome is displayed on Fig. [21 
we still compare the same quantities. Absolute errors on wave functions (left side) 
are much bigger for in this case. In particular, no new frequencies appear in 




Figure 2. Same as Fig. [U but cr = 2 and a = 2. (Critical case) 

the numerical solutions. The nonlinear effect boils down to a small change on the 
modulus of . 

4.3. Supercritical case. We close this first series of tests by considering a — 2 
and a = 1.5 as shown in Fig. [31 Of course, as no pointwise convergence is expected 
in this case, absolute errors are even bigger for both wave functions (left side) and 
moduli (right side). Of course, the size of the error on the modulus is much bigger 




Figure 3. Same as Fig. [1] but a — 2 and a = 1.5. (Supercritical case) 

too, and one should be extremely careful about the credit to give to the numerical 
simulations in the supercritical case. Indeed, this is a regime where a small error 
can be amplified at leading order (see [TTl ll2j V 

5. Visualization of the scattering operator 

We aim now at illustrating the results on scattering theory through numerical 
computations still achieved through time-splitting FFT schemes. The algorithm 
we used for the approximation of the scattering operator is based on a nonlinear 
time-splitting routine flanked by two free evolution steps (implemented the way 
recalled in the previous section): 

^V-f/(-r)oC/wL(2T)oC/(-T)V, r>i, 
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with U, Unl standing for the solution operators of equations ()2.ip and ()2.6p in 1-D 
with £ = 1 respectively. We used T = 55 in the computations hereafter. 

As no small parameter e is present in the problem, one may think that no major 
obstacle exists in carrying out this program; this isn't correct as the free evolutions 
can (and do!) dramatically increase the size of the computational domain for large 
T. It is interesting to notice that, in case one wants to use FFT-based schemes, 
both the direct computation for small e and the scattering operator approximation 
lead to a "large computational domain difficulty" : in the Fourier space for the first 
case, in the usual space for the second. 

A way to understand the scattering operator is to visualize the average effects 
of the nonlinearities appearing in equations of the form 

(5.1) idtux + ^dlux^ Xluxl^^ux, ux{t = 0,x) = cxp{-5x^), 

for various values of ct > 1 and A. Intuitively, as a increases, the nonlinearity 
becomes shorter range. Similarly, as A increases, the nonlinearity becomes stronger, 
and it should take a larger amount of time before we can consider it has become 
negligible. In all the tests we performed, it was somehow surprising to observe how 
fast the algorithm converges: one does not have to consider "very large" values of 
T so that 

Ui-T) o [/jvl(2T) o U{~T)i^ 
becomes stable and visually independent of T. 

5.1. Quintic nonlinearity (ct — 2). The parameter A controls in some sense the 
strength of the nonlinearitjQ inside (|5.ip . as can be seen on Fig. 01 This figure 
displays the position density of the initial data, the scattered solution for T = 
55 and a "mixed state" ux{t = 0) = U{-T)UNL{T)ux{t = 0). As our time- 
splitting/FFT algorithm preserves only the norm, but not the Hamiltonian, we 
first restricted ourselves to moderate values of A > (defocusing case). However, 




Figure 4. Initial data, numerical solution at T = 0, scattering 
state for ((57T|) with ct = 2 and A = 1 (left), A = 5 (right) 

as a numerical experiment, we wanted to display the outcome of our scheme for 
the stronger case A = 25 on Fig. [5l notice the change of shape in the scattered 
solution. Moreover, on this figure, we also tried to show what happens for A = — 1, 
that is to say for the focusing case despite there may be finite time blow-up (but 
there is scattering for small data). We checked that the energy associated to this 



through the stiffness of the associated differential equation. 
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data is (and remains) positive, a case where the virial identity , [14j . does not imply 
blow-up. The computational domain for these runs was [— IOOtt, IOOtt] and 2^^ — 1 
Fourier modes were used. 




Figure 5. Same as Fig.H but A = 25 (left), A = -1 (right) 



5.2. Power 3 (ct = 1.5). Now let's observe the effects of lowering the a value 
while keeping other parameters equal, see Fig. [S] It is interesting to see that the 
change of shape appearing for A = 25 is stronger than in the preceding case. On the 
contrary, the increase of the numerical solution's support is slightly less important. 
This hints that increasing the a value tends to expand the support of the scattered 




Figure 6. Same as Fig.lH but a = 1.5 and A = 5 (left), A = 25 (right) 

solution whereas increasing the A > value (defocusing case) leads to an oscillatory 
behavior. However, we stress that since the energy, 

E{t) := \\\d,u\t)\\l. + -h-\\u\t)f^+^, = E{Q), 

of the numerical solution changes more with a bigger A (its mass being always kept 
constant), these oscillations might be spurious. We actually don't know how this 
fact can be decided; our profiles have been checked to be stable on a finer grid. 

5.3. Power 6 (cr = 3). In order to get some numerical evidence about the de- 
pendence of the scattered solution on a, we display on Fig. [7] the outcome for 
(7 = 3. It is quite clear that the scattered solutions for both values of A are less 
peaked. Their support is bigger and the oscillations for A = 25 are weaker, their 
frequency remained the same though. This agrees with the behavior we sketched 
in the preceding subsection as a and A vary. 
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Figure 7. Same as Fig.gl but ct = 3 and A = 5 (left), A = 25 (right) 

6. Experiments on a cusp caustic 

Let us now go back to comparing the quadratic observables generated by numer- 
ical approximation of equations (|2.6p and (j2.ip endowed with a small parameter e 
in 1-D. In this section we fixed e ~ 1/150. Figure [5] displays the position density 
of the initial data for both equations, i.e. 

u'^{t = 0, x) = v^{t = 0,2:) = exp ^ - 2(a; - tt)^ - icos{x)/ej , x e [0, 27r], 

together with the position density of the numerical approximations of (|2.6p . (|2.ip 
in T = 3.5. The point here is to investigate what happens for the case of such a 
self-interfering Gaussian pulse, since no scattering theory is known for this problem. 
What we would like to check is whether the theoretical results on the focus point 
recalled and visualized in the preceding sections can be thought of as a guideline for 
this more complex case involving a non-trivial caustic. We shall observe position 
densities for the unique value of tr = 4 as a similar behavior has been seen to hold 
for different nonlinearities with convenient values of a. 4095 Fourier modes have 
been used in order to produce these results. 

6.1. Subcritical picture: a = 4. This case could be referred to as subcritical 
since it is noticeable on the top of Fig. [8] that the free and the nonlinear numerical 
solutions do agree for this reasonably small value of e. In particular, the frequen- 
cies of oscillations are identical. This is very similar compared to the behavior 
investigated in :8j. 

6.2. Critical picture: a — 3. The parameter a is now in a "critical range" as we 
observe that both solutions differ much more, but the frequency of the oscillations 
looks like being still the same in both cases. In order to establish this fact, we display 
on the left of Figl9]the FFT of the position densities: a peak at the same frequency 
is clearly noticeable. The nonlinear effect manifests itself through a change of order 
zero in the moduli, as we already observed on the right side of Fig. O notice also 
the similarity with the scattering state shown on Fig. [5] (right). This does agree 
with the ttc value derived in Section [2.21 

6.3. Supercritical picture: a = 2. In this last case, there is no similarity no 
more between the approximate solutions of (|2.6p . (|2.ip . as seen on both Fig. [8] and 
[9l Especially, the right side of Fig. [9] reveals that new frequencies show up inside 
the position density of the nonlinear solution. We have therefore a change of order 
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Figure 8. Position densities in the cusp caustic: 
to bottom) 



a = 4, 3, 2 (top 
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Figure 9. Fourier transforms of the positfon densities with a — 3 
(left), a = 2 (right) 



zero in the moduli and in the frequency. This is of course reminiscent of Fig. [3] in 
which a frequency doubling seems to show up in the supercritical regime. 

7. Conclusion 

We have presented the semi-classical limit for the nonlinear Schrodinger equation 
in the presence of a caustic. When the caustic is reduced to a point, the numerical 
experiments are in good agreement with the analytical results as far as the notion 
of criticality is concerned. However in the critical case, described by a nonlinear 
scattering operator, the leading order nonlinear effects are rather hard to visualize 
in the semi-classical limit. This is why we simulated the scattering operator in a 
separate way. 

Our numerical tests give encouraging evidence of new phenomena concerning the 
phase of the wave in the supercritical case when a focal point is formed (appearance 
of new frequencies). In the presence of a cusp caustic, the numerical experiments 
are in good agreement with the heuristic arguments that we presented here, for 
which no rigorous justification is available so far. 
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